!pip install \
numpy \
scipy \
statsmodels \
plotly \
numba
import numpy as np
import scipy as sp
import scipy.stats
import scipy.special
import scipy.integrate
import statsmodels as sm
from statsmodels.distributions.empirical_distribution import ECDF
import plotly
import plotly.graph_objs as go
from numba import njit
plotly.offline.init_notebook_mode()
Гипергеометрическое распределение - дискретное распределение, описывающее вероятность события, при котором ровно $k$ из $n$ случайно выбранных элементов окажутся помеченными, при этом выборка осуществляется из множества мощности $N$, в котором присутствует $m$ помеченных элементов. Считается, что каждый из элементов может быть выбран с одинаковой вероятностью $\frac{1}{N}$. Запишем это формально: $$\begin{gathered} N \in \mathbb{N},\ m \in \overline{0, N},\ n \in \overline{0, N},\\ k \in \overline{\max(0, m + n - N), \min(m, n)} \end{gathered}$$
Тогда $HG(N, m, n)$ описывает вероятность события, при котором ровно $k$ из $n$ элементов выборки окажутся помеченными: $$\left\{\ \xi \sim HG(N, m, n)\ \right\} \iff \left\{\mathbb{P}\left(\, \xi=k \,\right) = \frac{\binom{m}{k}\binom{N-m}{n-k}}{\binom{N}{n}}\right\}$$
По определению, математическое ожидание случайной величины – это ее $1^\text{й}$ начальный момент. Для начала, найдем $k^\text{й}$ начальный момент для $\xi$ (это понадобится для дальнейших выводов): $$\mathbb{E}\left[\, \xi^r \,\right] = \sum_{k=0}^{n} k^r \cdot \mathbb{P}\left(\, \xi=k \,\right) = \sum_{k=0}^{n} k^r\frac{\binom{m}{k}\binom{N-m}{n-k}}{\binom{N}{n}}$$ Можем считать, что сумма берется при $k=\overline{1,n}$, так как слагаемое при $k=0$ будет равно $0$. Заметим, что $$\begin{aligned} k\binom{m}{k} &= k \frac{m!}{k!(m-k)!} =\\ &= k \frac{m \cdot (m-1)!}{k \cdot (k-1)! \cdot (m-k)!} =\\ &= m \frac{(m-1)!}{(k-1)! \cdot \left(m-1 - (k-1)\right)!} =\\ &= m \binom{m-1}{k-1} \end{aligned}$$ и, как следствие, $$\binom{N}{n} = \frac{1}{n} \cdot n \cdot \binom{N}{n} = \frac{1}{n} N \binom{N-1}{n-1}$$ Подставим: $$\mathbb{E}\left[\, \xi^r \,\right] = \frac{n \cdot m}{N} \sum_{k=1}^{r-1} \frac{\binom{m-1}{k-1}\binom{N-m}{n-k}}{\binom{N-1}{n-1}}$$ Положим $j := k-1$ и изменим индекс суммирования на $j = \overline{0, n-1}$. Заметим, что $n - k = n - (j+1) = (n-1) - j$ и $N - m = (N-1) - (m-1)$: $$\mathbb{E}\left[\, \xi^r \,\right] = \frac{n \cdot m}{N} \color{dodgerblue}{\sum_{j=0}^{n-1} (j+1)^{r-1} \frac{\binom{m-1}{j}\binom{(N-1) - (m-1)}{(n-1) - j}}{\binom{N-1}{n-1}}}$$ Заметим, что выделенная часть выражения может быть записана, как $\mathbb{E}\left[\, (\theta+1)^{r-1} \,\right]$, где $\theta \sim HG(N-1, m-1, n-1)$. Следовательно, $$\mathbb{E}\left[\, \xi^r \,\right] = \frac{n \cdot m}{N} \mathbb{E}\left[\, (\theta+1)^{r-1} \,\right]$$ Таким образом, $$\boxed{ \mathbb{E}\left[\, \xi \,\right] = \frac{n \cdot m}{N} }$$
По определению дисперсии, $$\mathrm{Var}\left[\, \xi \,\right] = \mathbb{E}\left[\, \left(\ \xi - \mathbb{E}\left[\, \xi \,\right]\ \right)^2 \,\right] = \mathbb{E}\left[\, \xi^2 \,\right] - \left(\ \mathbb{E}\left[\, \xi \,\right]\ \right)^2$$ Выведем $2^\text{й}$ начальный момент: $$\mathbb{E}\left[\, \xi^2 \,\right] = \frac{n \cdot m}{N}\mathbb{E}\left[\, \theta+1 \,\right] = \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1}+1\right)$$ Подставим: $$\begin{aligned} \mathrm{Var}\left[\, \xi \,\right] &= \mathbb{E}\left[\, \xi^2 \,\right] - \left(\mathbb{E}\left[\, \xi \,\right]\right)^2 =\\ &= \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1}+1\right) - \left(\frac{n \cdot m}{N}\right)^2=\\ &= \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1} + 1 - \frac{n \cdot m}{N}\right) \end{aligned}$$ Таким образом, $$\boxed{ \mathrm{Var}\left[\, \xi \,\right] = \frac{n \cdot m}{N}\left(\frac{(n-1)(m-1)}{N-1} + 1 - \frac{n \cdot m}{N}\right) }$$
По определению, производящая функция вероятностей $G(z, \xi)$ – это математическое ожидание новой случайной величины $z^\xi$. То есть: $$G_\xi(z) = \mathbb{E}\left[\, z^\xi \,\right]$$ Для $\xi \sim HG(N, m, n)$ производящая функция выглядит так: $$G_\xi(z) = \binom{N-m}{n}\left({}_{2}F_{1}(-m, -n; N-m-n+1; z) - 1\right)$$ Здесь ${}_{2}F_{1}$ - это гипергеометрическая функция, определенная следующим образом: $${}_{2}F_{1}(a,b;c;z) = \sum_{n=0}^{\infty} \frac{a^{(n)} b^{(n)}}{c^{(n)}} \frac{z^n}{n!}$$ , а $x^{(n)}$ - возрастающий факториал, определенный как: $$x^{(n)} = \prod_{k=0}^{n-1} (x + k)$$
По определению, характеристическая функция случайно величины $\xi$ задается следующим образом: $$\varphi_\xi (t) = \mathbb{E}\left[\, e^{it\xi} \,\right]$$ Для $\xi \sim HG(N, m, n)$ характеристическая функция выглядит так: $$M_\xi(t) = \frac{\binom{N-D}{n}}{\binom{N}{n}}\ {}_{2}F_{1}\left(-n, -D; N-D-n+1; e^{it}\right)$$ Здесь ${}_{2}F_{1}$ - это гипергеометрическая функция.
Гистограмма – это графическое представление функции, приближающей плотность вероятности распределения на основе выборки из него.
Чтобы построить гистограмму, сначала нужно разбить множество значений выборки на несколько отрезков. Чаще всего, берут отрезки одинаковой длины, чтобы облегчить восприятие получившегося результата, однако это необязательно. Далее подсчитывается количество вхождений элементов выборки в каждый из отрезков и рисуются прямоугольники, по площади пропорциональные количеству попавших элементов выборки в соответствующий отрезок.
Вообще говоря, гистограмму можно использовать не только для приближения плотности на основе выборки, но и для визуализации самой плотности распределения, зная его плотность.
Мы будем строить гистограмму вероятностей, писать будем на языке Python3 и использовать следующие библиотеки:
Итак, для начала, определим класс гипергеометрического распределения
HG, который будет содержать в себе информацию о параметрах $N$, $m$ и
$n$ и предоставлять метод p(k), возвращающий вероятность принятия
случайной величиной значения k при данных параметрах:
class HG(object):
def __init__(self, N: int, m: int, n: int):
self.N = N
self.m = m
self.n = n
def p(self, k: int) -> float:
return sp.special.comb(self.m, k) * sp.special.comb(self.N-self.m, self.n-k) / sp.special.comb(self.N, self.n)
@property
def domain(self):
return (max(0, self.m + self.n - self.N), min(self.m, self.n))
def __str__(self) -> str:
return f'HG({self.N}, {self.m}, {self.n})'
Далее создадим объект случайной величины $\xi \sim HG(30, 15, 20)$:
xi = HG(80, 15, 40)
Следующим шагом, определим интервал $\overline{0, n}$, на котором мы будем рисовать нашу гистограмму:
hg_data_x = np.arange(*xi.domain)
И, наконец, построим гистограмму и выведем ее:
theoretical_hg_PDF = np.apply_along_axis(xi.p, 0, hg_data_x,)
theoretical_hg_PDF_trace=go.Scatter(
name='PDF',
legendgroup='PDF',
x=hg_data_x,
y=theoretical_hg_PDF,
mode='markers',
marker_color='blue',
)
hg_hist_fig = go.Figure(
data=(theoretical_hg_PDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\xi \sim ' + str(xi) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$\mathbb{P}(\xi=k)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$k$',
),
),
),
)
plotly.offline.iplot(hg_hist_fig)
По определению, функция распределения $F_\xi(k) = \mathbb{P}\left(\, \xi < k \,\right)$. Для дискретной случайной величины событие $\{\xi < k\} = \bigcup\limits_{i=0}^{k-1}\{\xi=i\}$. Каждое из событий $\{\xi=i\}\; \forall i \in \overline{0, k-1}$ являются попарно несовместными. То есть $\forall i,j \in \overline{0, k-1}: i \neq j$ выполняется $\{\xi=i\}\cap\{\xi=j\}=\emptyset$. Из этого следует, что $$\mathbb{P}\left(\, \xi < k \,\right) = \sum_{i=0}^{k-1}\mathbb{P}\left(\, \xi = i \,\right)$$ Подставим и получим: $$F_\xi(k) = \sum_{i=0}^{k-1}\mathbb{P}\left(\, \xi = i \,\right) = \sum_{i=0}^{k-1}\frac{\binom{m}{i}\binom{N-m}{n-i}}{\binom{N}{n}}$$
Построим график этой функции, учитывая, что аргументом $k$ должно быть натуральное число, не превосходящее $n$:
theoretical_hg_CDF = np.cumsum(np.apply_along_axis(xi.p, 0, hg_data_x))
theoretical_hg_CDF_trace=go.Scatter(
name='CDF',
legendgroup='CDF',
x=hg_data_x,
y=theoretical_hg_CDF,
line=go.scatter.Line(
shape='hv',
color='blue',
),
)
hg_dist_fig = go.Figure(
data=(theoretical_hg_CDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\xi \sim ' + str(xi) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$\mathbb{P}(\xi<k)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$k$',
),
),
),
)
hg_dist_fig.show()
Нормальное распределение - непрерывное распределение, описывающее поведение величины отклонения измеряемого значения $x$ от истинного значения $\mu$ (которое является математическим ожиданием) и в рамках некоторого разброса $\sigma$ (среднеквадратичного отклонения). Запишем это формально: $$\left\{ \eta \sim N(\mu, \sigma^2) \right\} \iff \left\{\begin{gathered} F_\eta(x) = \mathbb{P}\left(\, \eta < x \,\right) = \int_{-\infty}^{x} f_\eta(x)dx,\\ \text{где } f_\eta(x) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} \text{– плотность вероятности} \end{gathered}\right\}$$
Найдем математическое ожидание $\eta \sim N(\mu, \sigma^2)$: $$\begin{aligned} \mathbb{E}\left[\, \eta \,\right] &= \int_{-\infty}^{+\infty} x \cdot f_\eta(x)dx =\\ &= \int_{-\infty}^{+\infty} xe^{-\frac{(x-\mu)^2}{2\sigma^2}}dx =\\ &= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{+\infty} xe^{-\frac{(x-\mu)^2}{2\sigma^2}}dx \end{aligned}$$ Сделаем замену $t = \frac{x-\mu}{\sqrt{2}\sigma}$: $$\begin{aligned} \mathbb{E}\left[\, \eta \,\right] &= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{+\infty}(\sigma\sqrt{2}t + \mu) e^{-t^2} d\left(\frac{x-\mu}{\sqrt{2}\sigma}\right) =\\ &= \frac{\sigma\sqrt{2}}{\sqrt{\pi}}\int_{-\infty}^{+\infty}te^{-t^2}dt + \frac{\mu}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-t^2}dt =\\ &= \frac{\sigma\sqrt{2}}{\sqrt{\pi}}\left(\int_{-\infty}^{0}te^{-t^2}dt - \int_{-\infty}^{0}te^{-t^2}dt\right) + \frac{\mu}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-t^2}dt =\\ &= \frac{\mu}{\sqrt{\pi}}\int_{-\infty}^{+\infty}e^{-t^2}dt \end{aligned}$$ Заметим, что получившееся выражение содержит интеграл, который может быть сведен к интегралу Эйлера-Пуассона: $$\int_{-\infty}^{+\infty}e^{-t^2}dt = 2\int_{0}^{+\infty}e^{-t^2}dt = \sqrt{\pi}$$ Таким образом, $$\boxed{ \mathbb{E}\left[\, \eta \,\right] = \mu }$$
Сделаем ту же замену переменной $t = \frac{x-\mu}{\sqrt{2}\sigma}$, тогда $x = t\sqrt{2}\sigma+\mu$ и: $$\begin{aligned} \mathrm{Var}\left[\, \eta \,\right] &= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{+\infty}(\sqrt{2}\sigma)^2 t^2 e^{-t^2}d(t\sqrt{2}\sigma+\mu) =\\ &= \frac{2\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}t^2 e^{-t^2}dt \end{aligned}$$ Проинтегрируем по частям: $$\begin{aligned} \mathrm{Var}\left[\, \eta \,\right] &= \frac{\sigma^2}{\sqrt{\pi}}\int_{-\infty}^{+\infty}t 2t e^{-t^2} dt =\\ &= \frac{\sigma^2}{\sqrt{\pi}}\left(\left. -t e^{-t^2} \right|_{-\infty}^{+\infty} + \int_{-\infty}^{+\infty}e^{-t^2}dt\right) \end{aligned}$$ Здесь снова появляется интеграл Эйлера-Пуассона и, в итоге, получаем: $$\boxed{ \mathrm{Var}\left[\, \eta \,\right] = \sigma^2 }$$ То есть, $\sigma$ является среднеквадратичным отклонением.
Характеристическая функция для $\eta \sim N(\mu, \sigma^2)$ имеет вид: $$\varphi_\eta(t) = \exp\left(\mu it + \frac{\sigma^2 t^2}{2}\right)$$
Определим класс нормального распределения N, который будет содержать в
себе информацию о параметрах $\mu$ и $\sigma$ и предоставлять следующие
методы:
f(x) - возвращает значение плотности в точке xp(k) - возвращает
$\mathbb{P}\left(\, \eta < x \,\right) = \int_{-\infty}^x f_\eta(x)dx$class Norm(object):
def __init__(self, mu: float, sigma: float):
self.mu = mu
self.sigma = sigma
def f(self, x: float) -> float:
return np.exp(-((x-self.mu)**2)/(2*self.sigma**2))/(self.sigma*(2*np.pi)**.5)
def p(self, x: float) -> float:
return sp.integrate.quad(self.f, -np.inf, x)[0]
def __str__(self):
return f'N({self.mu}, {self.sigma}^2)'
Далее создадим объект случайной величины $\xi \sim HG(30, 15, 20)$:
eta = Norm(0, 1)
Следующим шагом, руководствуясь правилом трех сигм, определим интервал $(-3\sigma, 3\sigma)$, в котором окажутся все значения случайной величины с вероятностью более $0.99$:
norm_data_x = np.linspace(-3*eta.sigma, 3*eta.sigma, 100)
И, наконец, построим плотность, используя метод N.f нашего класса:
theoretical_norm_PDF = np.vectorize(eta.f)(norm_data_x)
theoretical_norm_PDF_trace = go.Scatter(
name='PDF',
legendgroup='PDF',
x=norm_data_x,
y=theoretical_norm_PDF,
line=go.scatter.Line(
color='blue',
),
)
norm_dens_fig = go.Figure(
data=(theoretical_norm_PDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\eta \sim ' + str(eta) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$f_\eta(x)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$x$',
),
),
),
)
plotly.offline.iplot(norm_dens_fig)
По определению, функция распределения $F_\eta(x) = \mathbb{P}\left(\, \eta < x \,\right)$. Для непрерывной случайной она определяется как интеграл от функции плотности вероятности: $$\mathbb{P}\left(\, \eta < x \,\right) = \int_{-\infty}^x f_\eta (x) dx$$. Подставим и получим: $$F_\xi(k) = \int_{-\infty}^x \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} dx = \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^x e^{-\frac{(x-\mu)^2}{2\sigma^2}} dx$$
Построим график этой функции, используя метод Norm.p:
theoretical_norm_CDF = np.vectorize(eta.p)(norm_data_x)
theoretical_norm_CDF_trace = go.Scatter(
name='CDF',
legendgroup='CDF',
x=norm_data_x,
y=theoretical_norm_CDF,
line=go.scatter.Line(
color='blue',
),
)
norm_dist_fig = go.Figure(
data=(theoretical_norm_CDF_trace,),
layout=go.Layout(
title=go.layout.Title(
text=r'$\eta \sim ' + str(eta) + '$',
x=.5,
),
yaxis=go.layout.YAxis(
title=go.layout.yaxis.Title(
text=r'$\mathbb{P}(\eta<x)$',
),
),
xaxis=go.layout.XAxis(
title=go.layout.xaxis.Title(
text=r'$x$',
),
),
),
)
plotly.offline.iplot(norm_dist_fig)
Основные алгоритмы моделирования брались из книги Р. Н. Вадзинского "Справочник по вероятностным распределениям"
Объявим константы:
QS = (5, 10, 10**2, 10**3, 10**5)
N = 5
Добавим к нашему классу HG метод gen, который будет возвращать реализацию из $HG(N, m, n)$:
@njit
def gen_hg(N: int, m: int, n: int) -> int:
x = 0
for _ in range(n):
if np.random.rand() < m/N:
x += 1
m -= 1
N -= 1
return x
HG.gen = lambda self: gen_hg(self.N, self.m, self.n)
Поясню происходящее: x – это количество помеченных элементов из тех, что мы выбрали.
Изначально x = 0, потому что мы еще не выбрали ни одного элемента.
Далее начинаем выбирать n элементов. Каждый раз, вероятность того, что мы выберем помеченный элемент равна m/N. Мы используем функцию numpy.random.rand, которая возвращает действительное число в интервале $[0, 1)$, и проверяем, меньше ли значение, чем m/N. Если это значение меньше, это означает, что мы выбрали один из помеченных элементов. В этом случае, мы увеличиваем x на единицу и уменьшаем количество оставшихся помеченных элементов на 1, потому что выборка ведется без возвращения.
Также, вне зависимости от того, был ли выбранный объект помеченным, мы уменьшаем количество объектов, из которых ведется выборка, на единицу.
Этот процесс повторяется ровно n раз, так как всего должно быть выбрано n элементов.
Для ускорения работы функции (она будет часто вызываться при построении выборок), мы используем декоратор из прекрасной библиотеки numba, который при первом вызове функции скомпилирует нашу функцию в машинный код, и, при последующих вызовах, будет вызываться уже скомпилированная функция, которая выполняется во много раз быстрее (в нашем случае прирост в скорости будет около 2000%).
Объем выборки будем обозначать $q$, чтобы не путать с параметром распреления $n$.
Напишем новый метод HG.sample_hg(q), который будет возвращать выборку объема $q$ из $HG(N, m, n)$:
def sample_hg(N: int, m: int, n: int, q: int) -> np.ndarray:
return np.fromiter((gen_hg(N, m, n) for _ in range(q)), np.int)
HG.sample = lambda self, q: sample_hg(self.N, self.m, self.n, q)
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ сгенерируем по 5 выборок $X_q^i$ объема $q$, где $i\in\overline{1, 5}$
samples_hg = {}
for q in QS:
samples_hg[q] = np.ndarray(shape=(N, q), dtype=np.int)
for i in range(N):
samples_hg[q][i] = xi.sample(q)
Теперь samples_norm – словарь, где ключом является объемом выборки, а значением по соответствующему ключу – массив из N = 5 выборок соответствующего объема из распределения $HG(80, 15, 20)$.
Выведем выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(samples_hg[q]):
print(f'#{i+1}: {s}')
Для каждой выбоки вычислим эмпирическую функцию распределения (ECDF) с помощью функции statsmodels.distributions.empirical_distribution.ECDF:
HG_ECDFs = {}
for q in samples_hg:
HG_ECDFs[q] = np.ndarray(shape=(samples_hg[q].shape[0], hg_data_x.shape[0]), dtype=np.float)
for i, sample in enumerate(samples_hg[q]):
HG_ECDFs[q][i] = ECDF(sample)(hg_data_x)
Теперь для каждой выборки построим эмпирическую функцию распределения, а также сравним получившуюся функцию с теоретической функцией распределения:
HG_ECDFs_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes='all',
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$k$',
y_title=r'$\mathbb{P}\left(\,\xi\leq k\,\right)$')
for i, q in enumerate(HG_ECDFs):
for j, sample in enumerate(HG_ECDFs[q]):
theoretical_hg_CDF_trace.showlegend = (i + j == 0)
HG_ECDFs_fig.add_trace(theoretical_hg_CDF_trace, row=i+1, col=j+1)
HG_ECDFs_fig.add_trace(
go.Scatter(
name='ECDF',
x=hg_data_x,
y=HG_ECDFs[q][j],
line=go.scatter.Line(
shape='hv',
color='red',
),
legendgroup='CDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
HG_ECDFs_fig.update_layout(height=800, width=900)
HG_ECDFs_fig.show()
Все в порядке, ECDF стремится к CDF при увеличении объема выборки.
Для каждого $q\in\{5, 10, 100,1000, 10^5\}$ найдем верхнюю границу разности каждой пары эмпирических функций распределения $\Delta_q^{i,j}$: $$\Delta_q^{i,j} = \sup\limits_{y\,\in\,\rm{supp}\,\xi}\left(\,\hat{F}\left(y, \overrightarrow{x_i}\right) - \hat{F}\left(y, \overrightarrow{x_j}\right)\,\right), \quad i,j\,\in\overline{0, 5}$$
for q in HG_ECDFs:
print(f'====== q = {q} ======:')
for i, ECDF_i in enumerate(HG_ECDFs[q]):
for j, ECDF_j in enumerate(HG_ECDFs[q][i+1:], start=i+1):
print(f'i = {i}; j = {j}; max delta = {np.absolute(ECDF_i - ECDF_j).max()}')
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим вариационный ряд выборки:
HG_VAR_ROWS = {}
for q in samples_hg:
HG_VAR_ROWS[q] = np.ndarray(shape=samples_hg[q].shape, dtype=np.int)
for i, sample in enumerate(samples_hg[q]):
HG_VAR_ROWS[q][i] = np.sort(sample)
Выведем получившиеся упорядоченные выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(HG_VAR_ROWS[q]):
print(f'#{i+1}: {s}')
Квантилью уровня $\alpha\in(0, 1)$ случайной величины $\xi$ называется такое число $x_\alpha\in\mathbb{R}$, что $$\begin{aligned} \mathbb{P}(\,\xi \leq x_\alpha\,) &\geq \alpha \\ \mathbb{P}(\,\xi \geq x_\alpha\,) &\geq 1-\alpha \\ \end{aligned}$$
Найдем квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ для гипергеометрического распределения $HG(80, 15, 40)$:
hg_quantiles = {}
for a in (.1, .5, .7):
for i, y in enumerate(theoretical_hg_CDF):
if y >= a:
hg_quantiles[a] = hg_data_x[i]
print(f'a = {a}, quantile = {hg_quantiles[a]}')
break
Выборочная квантиль уровня $\alpha\in(0, 1)$ выборки $\overrightarrow{X}$ - элемент вариационного ряда этой выборки стоящий на позиции $\left[\alpha\left|\overrightarrow{X}\right| + 1\right]$.
Для каждой выборки найдем квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ и сравним их с соответствующими квантилями данного распределения:
for q in QS:
print(f'======= q = {q} =======:')
for a in (.1, .5, .7):
print(f'------ a = {a} ------')
for i, sample in enumerate(HG_VAR_ROWS[q]):
quantile = sample[int(a*sample.shape[0])]
print(f'#{i}; quantile = {quantile:2}, real delta = {quantile - hg_quantiles[a]:2}')
Как и должно было быть, разность между выборочными квантилями и теоретическими стремится к нулю.
Вообще говоря, для вычисления квантилей выборок можно было воспользоваться функцией
np.quantile, но это было бы слишком просто...
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим гистограмму, а также сравним полученные графики с функцией вероятности $\mathbb{P}\left(\,\xi=k\,\right)$. Полигон частот строить не будем, потому что по сути, это просто альтернативный способ представления гистограммы, который будет нам только мешать анализировать графики.
HG_hists_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes=True,
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$k$',
y_title=r'$\mathbb{P}\left(\,\xi = k\,\right)$')
for i, q in enumerate(samples_hg):
for j, sample in enumerate(samples_hg[q]):
theoretical_hg_PDF_trace.showlegend = (i + j == 0)
HG_hists_fig.add_trace(theoretical_hg_PDF_trace, row=i+1, col=j+1)
HG_hists_fig.add_trace(
go.Histogram(
name='Histogram',
x=sample,
histnorm='probability',
marker_color='red',
legendgroup='PDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
HG_hists_fig.update_layout(height=800, width=900)
HG_hists_fig.show()
Как мы видим, закон больших чисел подтверждается графиками: с ростом объема выборки, значения гистограммы выборок все сильнее приближаются к истинному функции распределения.
С помощью формулы из справочника по вероятностным распределениям, мы можем получить сразу две реализации $x_i$ и $x_{i+1}$ из стандартного нормального распределения $N(0, 1)$, имея две реализации $r_i$ и $r_{i+1}$ из стандартного равномерного распределения на отрезке $[0, 1]$: $$\begin{aligned} x_i &= \sqrt{-2 \ln r_i} \sin(2\pi r_{i+1}) \\ x_{i+1} &= \sqrt{-2 \ln r_i} \cos(2\pi r_{i+1}) \end{aligned}$$
Чтобы получить реализацию $y_i$ из $N(\mu, \sigma)$: $$y_i = \mu + \sigma x_i$$
Для простоты, мы будем использовать только одну $y_i$ из двух возможных за раз:
@njit
def gen_norm(mu: float, sigma: float) -> float:
return mu + sigma * \
np.sqrt(-2 * np.log(np.random.rand())) * \
np.cos(2 * np.pi * np.random.rand())
Norm.gen = lambda self: gen_norm(self.mu, self.sigma)
Как говорилось ранее, объем выборки мы будем обозначать $q$.
Напишем метод HG.sample_norm(q), которsq будет возвращать выборку объема $q$ из $N(\mu, \sigma)$:
def sample_norm(mu: float, sigma: float, q: int) -> np.ndarray:
return np.fromiter((gen_norm(mu, sigma) for _ in range(q)), np.float)
Norm.sample = lambda self, q: sample_norm(self.mu, self.sigma, q)
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ сгенерируем по 5 выборок $X_\eta^{q,i}$ объема $q$, где $i\in\overline{1, 5}$
samples_norm = {}
for q in QS:
samples_norm[q] = np.ndarray(shape=(N, q), dtype=np.float64)
for i in range(N):
samples_norm[q][i] = eta.sample(q)
Теперь samples_norm – словарь, где ключом является объемом выборки, а значением по соответствующему ключу – массив из N = 5 выборок соответствующего объема.
Выведем выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(samples_norm[q]):
print(f'#{i+1}: {s}')
Для каждой выбоки вычислим эмпирическую функцию распределения (ECDF) с помощью функции statsmodels.distributions.empirical_distribution.ECDF:
Norm_ECDFs = {}
for q in samples_norm:
Norm_ECDFs[q] = np.ndarray(shape=(samples_norm[q].shape[0], norm_data_x.shape[0]), dtype=np.float)
for i, sample in enumerate(samples_norm[q]):
Norm_ECDFs[q][i] = ECDF(sample)(norm_data_x)
Теперь для каждой выборки построим эмпирическую функцию распределения, а также сравним получившуюся функцию с теоретической функцией распределения:
Norm_ECDFs_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes='all',
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$x$',
y_title=r'$\mathbb{P}\left(\,\eta\leq x\,\right)$')
for i, q in enumerate(Norm_ECDFs):
for j, sample in enumerate(Norm_ECDFs[q]):
theoretical_norm_CDF_trace.showlegend = (i + j == 0)
Norm_ECDFs_fig.add_trace(theoretical_norm_CDF_trace, row=i+1, col=j+1)
Norm_ECDFs_fig.add_trace(
go.Scatter(
name='ECDF',
x=norm_data_x,
y=Norm_ECDFs[q][j],
line=go.scatter.Line(
color='red',
),
legendgroup='CDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
Norm_ECDFs_fig.update_layout(height=800, width=900)
Norm_ECDFs_fig.show()
Тут опять все хорошо, значения ECDF стремятся к CDF при стремлении объема выборки к бесконечности.
Для каждого $q\in\{5, 10, 100,1000, 10^5\}$ найдем верхнюю границу разности каждой пары эмпирических функций распределения $\Delta_q^{i,j}$: $$\Delta_q^{i,j} = \sup\limits_{y\,\in\,\rm{supp}\,\eta}\left(\,\hat{F}\left(y, \overrightarrow{x_i}\right) - \hat{F}\left(y, \overrightarrow{x_j}\right)\,\right), \quad i,j\,\in\overline{0, 5}$$
for q in Norm_ECDFs:
print(f'====== q = {q} ======:')
for i, ECDF_i in enumerate(Norm_ECDFs[q]):
for j, ECDF_j in enumerate(Norm_ECDFs[q][i+1:], start=i+1):
print(f'i = {i}; j = {j}; max delta = {np.absolute(ECDF_i - ECDF_j).max()}')
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим вариационный ряд выборки:
Norm_VAR_ROWS = {}
for q in samples_norm:
Norm_VAR_ROWS[q] = np.ndarray(shape=samples_norm[q].shape, dtype=np.float64)
for i, sample in enumerate(samples_norm[q]):
Norm_VAR_ROWS[q][i] = np.sort(sample)
Выведем получившиеся упорядоченные выборки для q = 5 и q = 10:
for q in (5, 10):
print(f'====== q = {q} ======:')
for i, s in enumerate(Norm_VAR_ROWS[q]):
print(f'#{i+1}: {s}')
Найдем квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ для нормального распределения $N(0, 1^2)$:
norm_quantiles = {}
for a in (.1, .5, .7):
for i, y in enumerate(theoretical_norm_CDF):
if y >= a:
norm_quantiles[a] = norm_data_x[i]
print(f'a = {a}, quantile = {norm_quantiles[a]}')
break
Теперь для каждой выборки найдем выборочные квантили уровней $\alpha\in\{0.1, 0.5, 0.7\}$ и сравним их с соответствующими квантилями данного распределения:
for q in QS:
print(f'======= q = {q} =======:')
for a in (.1, .5, .7):
print(f'------ a = {a} ------')
for i, sample in enumerate(Norm_VAR_ROWS[q]):
quantile = sample[int(a*sample.shape[0])]
print(f'#{i}; quantile = {quantile:2}, real delta = {quantile - norm_quantiles[a]:2}')
Как и в случае с гипергеометрическим распределением, разность между выборочными квантилями и теоретическими стремится к нулю, что есть правильно.
Для каждого $q\in\{5, 10, 100, 1000, 10^5\}$ построим гистограмму, а также сравним полученные графики с плотностью распределения $f_\eta(x)$. Как и в соответствующем пункте про гипергеометрическое распределение, полигон частот строить не будем.
Norm_hists_fig = plotly.subplots.make_subplots(rows=len(QS),
cols=N,
shared_xaxes='all',
shared_yaxes=True,
row_titles=[f'q = {q}' for q in QS],
column_titles=[f'#{i}' for i in range(N)],
x_title=r'$x$',
y_title=r'$f_\eta(x)$')
for i, q in enumerate(samples_norm):
for j, sample in enumerate(samples_norm[q]):
theoretical_norm_PDF_trace.showlegend = (i + j == 0)
Norm_hists_fig.add_trace(theoretical_norm_PDF_trace, row=i+1, col=j+1)
Norm_hists_fig.add_trace(
go.Histogram(
name='Histogram',
x=sample,
histnorm='probability density',
marker_color='red',
legendgroup='PDF',
showlegend=(i + j == 0),
),
row=i+1,
col=j+1)
Norm_hists_fig.update_layout(height=800, width=900)
Norm_hists_fig.show()